Draft version May 21, 2009 

Preprint typeset using 1^1^^ style cmulatcapj v. 08/13/06 



O 

o 

(N 

>>: 



(N 

o 

u 

6 



> 

o 

rn 

vn 
o 
o\ 
o 



ANISOTROPIC AGN OUTFLOWS AND ENRICHMENT OF THE INTERGALACTIC MEDIUM 
Joel Germain^--, Paramita Barai^ ''^, and Hugo Martel^-^ 

Draft version May 21, 2009 

ABSTRACT 

We investigate the cosmological-scale influence of outflows driven by AGNs on metal enrichment 
of the intergalactic medium. AGNs are located in dense cosmological structures which tend to be 
anisotropic. We designed a semi-analytical model for anisotropic AGN outflows which expand away 
along the direction of least resistance. This model was implemented into a cosmological numeri- 
cal simulation algorithm for simulating the growth of large-scale structure in the universe. Using 
this modified algorithm, we perform a series of 9 simulations inside cosmological volumes of size 
(128 /i^^Mpc)'^, in a concordance ACDM universe, varying the opening angle of the outflows, the 
lifetimes of the AGNs, their kinetic fractions, and their level of clustering. For each simulation, we 
compute the volume fraction of the IGM enriched in metals by the outflows. The resulting enriched 
volume fractions are relatively small at z > 2.5, and then grow rapidly afterward up to z = 0. We find 
that AGN outflows enrich from 65% to 100% of the entire universe at the present epoch, for different 
values of the model parameters. The enriched volume fraction depends weakly on the opening angle 
of the outflows. However, increasingly anisotropic outflows preferentially enrich undcrdcnse regions, 
a trend found more prominent at higher redshifts and decreasing at lower redshifts. The enriched 
volume fraction increases with increasing kinetic fraction and decreasing AGN lifetime and level of 
clustering. 

Subject headings: cosmology: theory — galaxies: active — intergalactic medium — quasars: general 
— methods: N-body simulations. 



1. INTRODUCTION 

Active galaxies are believed to be powered by ac- 
cretion of matter onto the central supermassive black 
holes (SMBHs) (e.g., 'K ormendv fc Richstond Il995t 
iFerrarese fc FordI [2005), liberating enormous amounts 
of energy (often in the forms of ejected energetic 
jets/outflows), affecting their environments from pc to 
Mpc scales. Active Galactic Nuclei (AGN) are believed 
to influence the formation and evolution of galaxies and 
large-scale structures over the cosmic epoch, in the form 
of feedback, whereby the overall properties of a galaxy 
can be reg ulated by its central BH (e . g.. iSilk fc Reea 
1998; Ki ngI [200l fWvithe fc Loebl [200l iGranato et all 
20041: iMur rav et all 120051: iBege lman fc NathI 120051 : 
Pipino et a l. 2009; Ciotti et al. 2009|). Recent concor- 



dance models of galaxy formation through hierarchical 
clustering in the cold dark matter cosmology invoke feed- 
back from AGN to explain several observations (jBestI 
l2007l for a review), such as the central SMBH - host 
galaxy bulge correc t ions (e.g., Magorrian e t al.l 119981 : 
iGebhardt et all 120001: iMcLure fc Dunlop 2002yand the 
sharp cutoff at the bright end of the galaxy luminosity 
function. 

During the quasar era (1 < z < 3), the star- formation 
rate, the comoving density of AGN, and the merger rate 
of galaxies, are all observed to reach their peak values. 
Such a common trend points to a coevolution scenario: 
galaxies and AGN are argued to have evolved t o gether 
influencing each other's growth. iPolletta et al.l (|2008l ) 
discovered two sources at z ~ 3.5 exhibiting both pow- 
erful starburst and AGN activities, which is interpreted 
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as a coevolution phase of star formation and AGN as 
predicted by various models of galaxy formation and 
evolution. Ob servations of an AGN-starburst galaxy by 
ICroston et al.l ([20081 show that AGN feedback has a sub- 
stantial impact on the g alaxy and its surroundin gs. Cos- 
mological simulations bv lDi Matteo et all (|2005| ) indicate 
that the energy released from an accreting SMBH can 
halt further accretion onto the BH and drive away gas, 
self-regulating galaxy growth by shutting off the quasar 
and quenching star-formation in the galaxy, and thus 
reddening it rapidly. 

The nature of feedback from AGN on star/galaxy for- 
mation can be either negative (quenching) or positive 
(enhancement), as shown by different studies. The exact 
role played by AGNs possibly depends on one or more 
factors, such as radio-loudness, jet power, lifetime/duty 
cycle, environmental factors like ambient gas density, etc. 
In galaxy clusters, AGN outffows are believed to stop 
the cooling flow and heat up the intracluster medium by 
giant cavities, buoyant bu bbles and shock fronts (e.g., 
iMcNamara fc NulsenI [20071 ). At the same time, substan- 
tial star-formation rates are observed in the central radio 
galaxie s of cool ing flow clusters triggered by the radio 
source (jMcNama ra 2002i). 

Some observations (e.g., ISchawinski et al.l I2007D in- 
dicate that AGN quench star-formation in their host 
galaxies at z < 2 turning them red and dead, in the 
process generating the observed bimodal color distribu- 
tion of galaxies. Hydrodynamical simulations of ma- 
jor mergers of equal- and unequal-mass disk and ellip- 
tical galaxi es show that BH feedback can terminat e star- 
formation (ISpringel et al.l l2005l: [Johansson et al.l [2008) . 
lAntonuccio-Delogu fc Silkl ()2008l ) analyzed the different 
physical factors impacting the star-formation rate, and 
found that there is suppression of star-formation in 
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galaxies from AGN jet-induced feedback. Sem i-analytic 
models (|Bower et albooeHCroton et alJl2006D show that 
feedback due to AGN can quench cooling flows and 
star formation in massive halos, explaining the observed 
downsizing in the stellar populations of galaxies (e.g., 
lJuneau et aTll2005f ). 

Shock-induced and jet-induced star formation in radio 
galaxies can explain the radio-optical alignment obser- 
vations. In hydrodynamic simulations including radia- 
tive cooling, AGN-driven shock waves compress a gas 
cloud, causing it to break up into multiple smaller, dense 
fragments, which survive for many dynam ical timescales, 
turn Jeans unstable, and form stars (e.g..lMellema et al.l 
I2002t [Fragile et alll2004 Ivan Breugel et al.ll2004D . Ex- 
panding radio jets/lobes generate shocks which propa- 
gate through inhomogeneous clumpy ambient gas clouds 
and trigger gravitational collapse of overdense ambient 
gas clouds, and star f ormation (e.g.. 'Be gelman fc Cioffil 
[l98l IDe Yo{m3[l989l : IRccs 1989; Bickncl leral.ll2000n . 

Studies claim to have found direct observational evi- 
dence for AGN feedback aX, z ^ 2. iNesvadba fc Lehnert 
using VLT data, identified kpc-sized outflows of 
ionized gas in z ~ 2 — 3 radio galaxies. These bipolar 
outflows have the expected signatures of being powerful 
AGN-driven winds, energetic enough to terminate star 
formation in the massive host galaxies. At the same time 
there have been few observ ations (e.g., iKrongold et al.l 
120071: ISchlesinger et al.ll2009D implying that quasar winds 
are unlikely to produce important evolutionary effects on 
their larger environments, which poses a challenge for the 
scenario of strong AGN feedback in galaxy evolution. 

A fraction of AGN are radio loud, and there are claims 
that the radio-loud ness depends on SMBH mass a nd Ed- 
dington ratio (e.g.. lLaoH[2000l : iRafter et al.ll2009[ ). The 
lifetimes of radio sources 10^ — 10® yrs) are observa- 
tionally inferred to be signiflcantly shorter than the ages 
of their host galaxies, suggesting that the hosts cycle be- 
tween radio-loud and radio-quiet phases. AGN feedback 
is a natural explanation: the hot ISM in a galaxy episod- 
ically cools to fuel the central AGN, which triggers radio 
activity, and then the em anated radio jets reheat the gas 
(|Ciotti fc Ostrikeil[2007l) . 

A large fraction of AGN are obs erved to host outflows , 
in a wide variety of forms (see ICrenshaw et al.l l2003t 
lBegelmanll2004l: lEverettll2007l for reviews). Radio galax- 
ies with collimated relativistic jets and/or huge over- 
pressured cocoons emitting radio synchrotron r adiation 
constitute ~ 10% of all quasars (|Petersonlll997f ). Blue- 
shifted broad absorption lines (BALs) in the UV and op- 
tical are seen in an a dditional ~ 10— 15% of QSOs (e.g., 
iReichard et~ani2003f l. Many BAL and other quasars m 
SDSS exhibit highly-ioniz ed broad blue-s hifted emission 
lines of e.g., C IV (e.g., iRichards et"ani2002 ). Intrin- 
sic absorption lines in the UV have b een detected in 
> 50% of Seyfert I galaxies studied bv ICrenshaw et aH 
1X999). Outflows have been observed in [ QUI] emis- 
sion lines in nearby Seyfert galaxies (e.g., ^D as et aTl 
120051) . Warm absorbers seen in X-rays indicate ion- 
ized outflows in Seyferts a nd QSOs (e.g., iBlustin et al.l 
120051: IKrongold et al.|[2007l ). More high- velocity highly- 
ionized outflows have been detected as absorption lines 
in X-rays (e.g., 'Chartas et al.' 2003; Pou nds et aI1l2003l : 
iDasgupta et"al.i ,2005t ,0'Brien et al. 2005[k The full 
mechanism of AGN feedback and how it operates on dif- 



ferent scales is poorly understood, both in theoretical 
and observational aspects. In this paper we investigate 
the impact of energetic outflows emanating from AGN 
on large-scale cosmological volumes. 

There have been previous studies on the 

cosmological i mpact of quas ar outflows at 

large s cales (iFurlanetto fc Loebl l2001l hereafter 
FLOl; ISc annapieco fc__Oh| |2004j, here after SO04 ; 
iLevine fc G nedin 2005, hereafter LG05). iBarail (j2008l ) 
used cosmological simulations to investigate the large- 
scale influence of radio galaxies over the Hubble time. 
All these studies have considered outflows expanding 
with a spherical geometry. However, in realistic cosmo- 
logical scenarios, where the density distribution show 
significant structures in the form of filaments, pancakes, 
etc., outflows are expected to expand anisotropically on 
large scales. Recent observations support the picture 
of anisotropic b ipolar outflows on different scales. 
INesvadba et al.l (|2008f ) found spectroscopic evidence 
for bipolar outflows in three powerful radio galaxies at 
z ~ 2 — 3, with kinetic energies equivalent to 0.2% of 
the rest-mass of the SMBH. These outflows possibly 
indicate a significant phase in the evolution of the host 
galaxy . On small scales, lAnathpindika fc WhitworthI 
([2008') showed that the observed outfiows from young 
stellar objects tend to be orthogonal to the filaments 
that contain their driving s ources. Numerical s i mula- 
tions on galaxy-s cales (e.g., iMac Low fc Ferraral 119991 : 
iRecchi et al.ll2009l ) show that bipolar winds and outfiows 
expand preferentially along the steepest density slope, 
i.e., the direction perpendicular to the galaxy plane, 
while simulations at larger scales (|Martel fc Shapiro! 
l2001al |bi) show that outfiows emerging from cosmological 
structures expand preferentially along the direction of 
least resistance. 

The goal of our work is to study the global large- 
scale influence of outflows from AGN on the intergalactic 
medium (IGM) in a cosmological context. We have de- 
signed a semi-analytical model for anisotropic outflows, 
which we implemented into a numerical algorithm for 
cosmological simulations. Using this algorithm, we sim- 
ulate the propagation of AGN-driven outflows into the 
IGM, in a ACDM universe. We then explore the large- 
scale impact of the cosmological population of AGN out- 
flows over the age of the universe. In this paper, we 
focus on metal enrichment of the IGM, estimate the vol- 
ume fraction of the universe enriched by AGN outflows 
as a function of redshift, and its dependence on the var- 
ious parameters of the model. In a forthcoming paper, 
we will will focus on the metal content and metallicity of 
the IGM and its observational consequences. 

This paper is organized as follows. In fj^lwe describe 
our semi- analytical model for anisotropic outflows, and 
its implementation into a cosmological N-body simula- 
tion algorithm. The results are presented and discussed 
in fj3l We present our conclusions in ijH 

2. THE NUMERICAL METHOD 

Our numerical method consists of three ingredients, 
(1) a cosmological Particle- Mesh algorithm for simulat- 
ing the formation and evolution of large-scale structures 
in the universe, (2) a model for the masses, formation 
epochs, lifetimes, and spatial distributions of AGNs, 
and (3) a semi-analytical model for the propagation of 
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Fig. 1. — Redshift distributio n of a total of 1535 362 sources 
generated according to the QLF ( il2.2| l in various luminosity bins, 
with parameter values: Tagn = 100 Myr, and /outflow = 0.6. The 
number of sources A'^agn in each luminosity bin is indicated. 

anisotropic outflows and the resulting metal-enrichment 
of the IGM. We discuss these various ingredients in the 
following subsections. 

2.1. The PM Algorithm 

We simulate the growth of large-scale structure in 
a cubic cosmological volume of comoving size Lbox = 
128 Mpc = 182.6 Mpc with periodic bound- 

ary conditions, using a Par ticle-Mesh (PM) algorithm 
(jHocknev fc EastwoodI flQSl . We use 256'^ equal mass 
particles and a 512'^ grid. This corresponds to a par- 
ticle mass mpart = 1-38 x IO^^Mq, and a grid spacing 
A = 0.357 Mpc. Note that this length resolution is suf- 
ficient for our purpose; we do not need the extra length 
resolution that a P'^M algorithm would provide, and us- 
ing PM instead of P'^M results in a major speed-up of 
the calculation. LG05 also used a PM algorithm for their 
study of AGN outflows. 

We consider a ACDM model with a present baryon 
density parameter Slb^o = 0.0462, total matter (baryons 
+ dark matter) density parameter fip = 0.279, cosmo- 
logical constant 0.^^ = 0.721, Hubble constant iJp — 
70.1kms"iMpc"i [h = 0.701), primordial tilt = 
0.960, and CMB temperature Tcmb = 2.725, consis- 
tent with the results of WMAP5 combined with the data 
fro m baryonic acoustic o scillations and supernova stud- 
ies (jHinshaw et al.l[2008l ). We generate initial conditions 
at redshift z = 24, and evolve the cosmological volume 
up to a final redshift z = 0. 

2.2. Distribution in Redshift and Luminosity 

We determine the luminosities and birth times of the 
cosmological AGN population by adopting the rcdshift- 
dependent luminosity distribution of AGNs from the 



work of lHopkins et al.l (|2007f ) on bolometric quasar lumi- 
nosity function (QLF). The QLF is expressed as a stan- 
dard double power law. 



d<i> 



(1) 



dlogL {L/L^p^ +{L/L^p-' 

which gives the number of quasars per unit comoving 
volume, per unit log luminosity interval. The values of 
the amplitude 0*, break luminosity L^,, and the slopes, 
7i and 72 ev olve with redsh i ft. W e use the values given 
in Table 2 of lHopkins et al.l (|2007f ). which cover the red- 
shift range z = 6.0 — 0.1. For a given redshift z, we get 
the luminosity distribution by interpolating between two 
consecutive lines in that table. For redshifts z < 0.1, we 
simply use the luminosity function at z = 0.1. 

We assume t hat a fraction /outflow = 0-6 o f AGNs pro- 
duce outflows (jGangulv fc Brotherton|[2008f ). The num- 
ber of AGNs with outflows in the simulation box of co- 
moving volume Vbo 

terval [L.L + dL] at redshift z is given by 



box = .^box' ^^'^ within luminosity in- 



N{L,z) = /outfiowVbox'/'(i, z)(i[logL] 



(2) 



We assume that AGNs have luminosities between L,„in = 
lO^L© and L^ax = IQ^^Lq (Crenshaw et al. 2003). The 
fiducial value for the AGN activity lifetime is taken as 
Jagn = 10® yr, with other values considered in ii3.3l 

We calculate the number and birth time of the AGNs as 
follows. We divide the luminosity range L = [Lmin, ^max] 
in bins of size Alog(L/LQ) = 0.1. For each luminosity 
bin, we start at redshift z — Zmax — 6, and calculate the 
number of AGNs with outfiows in that bin at that red- 
shift using equation We then assign to each AGN a 
random age tagc between (AGN just being born) and 
2agn (AGN just about to die), and calculate the birth 
redshift Zbirth = z{tm&^ ~ ^agc), where tmax is the cos- 
mic time corresponding to redshift Zmax- We then move 
forward in time by intervals of At = 10 Myr.'^ For each 
new time t, we calculate the number of AGNs with out- 
fiows using equation ([2]), subtract from that number the 
number of AGNs present at earlier times that are still 
alive at that time, and assign to each new AGN an age 
between and At (since, if that AGN was not present at 
the previous time, it must have appeared during the last 
interval At). 

Using the QLF, we obtain the entire cosmological pop- 
ulation of AGNs in the simulation box starting from 
-Zmax, namely the birth redshift (zbir), switch-ofF redshift 
(zoff) and bolometric luminosity (iboi) of each source. 
Figure [1] shows the redshift distribution of the total pop- 
ulation of 1 535 362 sources produced using Tagn = 10® 

yr- 



Spatial Location of A GNs 
AGNs have been observed (e.g. 



2.3. 

By far the 

iD'Odorico et al.ll2008f ) to be hosted in high density re- 
gions of the universe. In order to determine the spatial 
location of the AGNs in our simulations, we consider 
the density distribution on the 512 x 512 x 512 grid, 
as calculated by the PM code. We filter this density 
using a ga ussian filter containing a m ass lO^^M^g (e.g., 
[Kauffmann et al.l 120031 : iHickox et al.l 12009.) , considering 

^ The length of this interval does not matter as long as it is 
sufficiently smaller than Tagn- 
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that as the minimum mass of a galaxy that might host 
an AGN (for det ails of the fil tering technique, we refer 
the reader to §5 of lMarte]||2005D . We then identify all grid 
points where the value of the filtered density exceeds the 
values at the 26 neighboring grid points. These arc the 
locations of the density peaks. 

At each timestep of the simulation, we spatially locate 
the new AGNs born during that epoch ( H2.21 whose Zbir 
values fall within the timestep interval) at the local den- 
sity peaks in the cosmological volume. We should have 
enough peaks to assign to each new AGN an unique lo- 
cation, but at the same time we should locate the AGNs 
preferentially in high-density regions. So we adopt a 
redshift-dependent limiting density pjjj^ to select poten- 
tial peaks for locating the sources. We consider the peaks 
that have a filtered density p > p^^^ , and each new AGN 
is located at the center of one such peak, selected ran- 
domly but excluding peaks already containing an AGN. 
The number of AGN increases with time, hence we need 
to reduce the value of piim with redshift in order to se- 
lect a sufficient number of peaks. We use p^^^ = lOp 
when z > 3, p^^^ = 7p when 3 > z > 2, p^^^^ = hp when 
2 > z > 1, Piijjj = 3p when 1 > z > 0.75, and py^^^ = 2p 
when z < 0.75, where p{z) = moH§{l + zf/SiiG is the 
mean density of the universe at that redshift. 

Note that our method for locating AGNs differs from 
the one used by LG05. In their simulations, AGNs could 
be located anywhere in the computational volume, with 
a probability that depended on the local density. We 
have to restrict the potential locations of AGNs to lo- 
cal density peaks, to be consistent with our anisotropic 
outflow model (see below). This is not a severe restric- 
tion, because we do expect massive galaxies that host 
AGNs to form at the locations of density peaks anyway. 
We consider alternate methods of locating the AGNs in 



2.4. Anisotropic Outflows 

iPieri et all (|2007f) (hereafter PMG07) have developed 
an anisotropic outflow model, in which outflows expand- 
ing into an anisotropic medium follows the path of least 
resistance. The model was applied to supernovae-driven 
outflows generated by low-mass galaxies. In this paper 
we apply the same model to AGN-driven outflows gen- 
erated by massive galaxies. The model is essentially the 
same, except for some additional terms in the equations 
driving the outflow. We refer the reader to PMG07 for 
details. 

The outflow is approximated as two expanding bipolar 
cones with opening angle a, expanding radially in oppo- 
site directions, as shown in Figure [H The opening angle 
is treated as a free parameter. The limits a = tt and 
a <C 1 correspond to the cases of isotropic outflows and 
jets, respectively. Note that in a given simulation, we use 
the same value of a for all outflows. The total volume of 
the outflow is 47ri?^(l - cosq;/2)/3. 

We assume that the outflow will follow the direction of 
least resistance as it travels away from the source. Our 
technique for d etermining thi s direction was presented 
in PMG07 and iGrenonI (|2007t ). Essentially, we perform 
a second-order Taylor expansion of the density around 




isotropic anisotropic 



Fig. 2. — Geometry of the outflow. R{t) and a are the radius 
and opening angle of the outflow, respectively. The unit vector 
e indicates the direction of least resistance. Figure taken from 
PMG07. 

each peak where a source is located: 

p{x,y, z) = ppcak-Ax^ -By^ -C z"^ -2Dxy-2Exz-2Fyz , 

(3) 

where the coordinates x, y, z are measured from the po- 
sition of the peak. Since the peak is by definition a local 
maximum, there are no linear terms in the expansion. In 
practice, we determine the coefficients A through F by 
performing a least-square fit of equation Q to all the 
grid points within a distance 2A from the peak (A being 
the grid spacing). We then rotate the coordinate axes 
such that, in the new coordinate system (x',y',z'), the 
cross-terms vanish: 

p{x', y', z') = ppeak - A'x'^ - B'y'^ - C z'^ . (4) 

The three coefficients A', B' , C are always positive, oth- 
erwise we would not have a peak, but rather a local min- 
imum or a saddle point. The largest of these coefficients 
gives us the direction along which the density drops the 
fastest as we move away from the peak. We take it as 
the direction of least resistance. 

To obtain a complete description of the outflow, we 
still need to provide an outflow model that will give 
us the time-dependent radius R{t). This is where our 
model differs from the one of PMG07, which is designed 
for supernovae-driven outflows. We describe our outflow 
model in the next subsection. 

2.5. Outflow Model 

Despite the observational differences between various 
kinds of outflows, the important point relevant for the 
present study is that the AGNs hosting outflows form 
a random subset of the whole AGN population (SO04). 
For simplicity we assume the same outflow model for all 
kinds of AGN, as in FLOl, SO04, and LG05. Neglect- 
ing outflow formation timescales (since we are bound by 
simulation resolutions), we consider that each AGN pro- 
duces an outflow right from its birth. The outflow ex- 
pands into the IGM with an anisotropic geometry, get- 
ting channeled into low-density regions. 

We set the kinetic luminosity Lk carried by the AGN 
jets equal to a constant fraction of the bolometric lumi- 
nosity, 

Lk — ex-^boi, (5) 

where tK is the kinetic fraction. F or our initial simu- 
lations, we use the value tK = 0.1 ([Willott et al.llT999l : 
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FLOl; LG05: IChartas et al.l[2007l : iShankar et~aI1[200l . 
which is probably an upper Hmit for BAL outflows 
()Nath fc Rovchowdhurvl I2002D . but we consider other 
values of ck in gMl SO04 used ex = 0.05. The total 
kinetic energy transported by the jets during an AGN's 
active lifetime, Ek = LkTag^, is assumed to be con- 
verted to the thermal energy of the outflow. 

We assume that the central black hole (of mass A/bh) 
radiates at the Eddington limit, and provides the AGN 
bolometric luminosity, 



AnGcm 



Eddington 



1.26 X lO^^ergs 



-1 



^Mbh 



Mbh 



(6) 



Then using the ratio of central BH mass to the 
galaxy bulge ste l lar m as s, MeHZ-^buige = 0.0 02 (e.g., 
Magorrian erall Il998t iGebhardt etall 120001 : FLOl; 



Marconi fc Hlmtl l2003'). and scaling it with the universal 



ratio of the density of matter to baryons, we obtain the 
total (baryonic+dark matter) mass of the galaxy hosting 
the AGN, Mgai = (f}o/f^b,o)A^buigo. 

The model of a spherical bubble expanding as a 
thin shell (of mass M s) in a cosmological volume (e.g., 
iTegmark et al.|[T993[ ) is used to obtain the radius, R, of 
the outflow. The rate at which IGM mass is swept up by 
an anisotropic outflow is given by 




cos (a/2)] 



R-Vr^ 



Vp> R; 
v„ < R: 



where px is the density of the external gas, and Vp is 
the gas velocity due to infall onto the dark matter halo 
hosting the AGN. As a simplifying approximation, we 
consider that the gas density is equal to the mean baryon 
density of the universe at the corresponding redshift. 



BvrG 



(l + z)= 



(8) 



For the gas infall velocity, we assume Vp = H{z)R, where 
H{z) is the Hubble constant at redshift z. This is simply 
the Hubble expansion of a cosmological volume. 

After correcting for anisotropic expansion, the acceler- 
ation of the shell can be written as, 

•• 47ri?2 [1 -cos(a/2)] , , 



G_ 



Md{R) + Afgal 



+nAiz)H^{z)R- 



r 

Ms 



R-Vp{R) 



(9) 



where px is the thermal pressure inside the outflow, pB is 
the magnetic pressure, p^ is the pressure of the external 
IGM, and Md{R) is the mass of matter lying inside the 
shell. In equation ([9]), the first term represents the pres- 
sure gradient driving the outflow outwards. The second 
term is the gravitational deceleration caused by matter 
existing inside the outflow radius (attraction between the 
outflow shell and the background halo + host galaxy), 
and the self gravity of the shell. The third term is an 
acceleration due to the cosmological constant. The final 



term is a drag force caused by sweeping up the IGM and 
accelerating it from velocity Vp to R. 

This formulation for the evolution of outflow radius 
[eq. ([5])] is almost same as that of FLOl, except that 
we have a factor of [1 — cos(a/2)] in the first term to 
take i nto account of t he an isotropic shape of the out- 
flow. ITegmark et all (|1993f) and PMG07 used a simi- 
lar equation, but they do not have the following terms 
which we include for AGN outflows: the magnetic pres- 
sure (ps), the deceleration due to self-gravity of the shell 
{GMs/2R^), and the acceleration due to the cosmologi- 
cal constant {Vlf^H^R). 

The external gas pressure is obtained using Px{z) = 
Px{z)kTx/ p, where the external temperature is fixed at 
Tx — 10^ K assuming a photoheated ambient medium, 
and p = 0.611 amu is the mean molecular mass, as- 
suming that the medium ahead of the outflow has been 
photoionized (PMG07). The working expression for the 
pressure of the IGM is then 



Px{z) 



'SHokTx 



zf 



(10) 



Making again the approximation that the density of 
the background through which the AGN outflows propa- 
gate is equal to the mean matter density of the universe, 
we obtain the halo mass enclosed within the shell, 

, 47ri?3 n(z)H^iz)R^ , , 

M4R, z) = -^P{z) = \q ■ (11) 

The thermal pressure is provided by the jet kinetic 
power when the AGN is active. It undergoes expansion 
losses and varies at a rate 



A 



Pt 



27ri?3 [1 - cos (a/2)] 



R 
R 



(12) 



The first and second terms represent the increase in pres- 
sure caused by injection of thermal energy into the out- 
flow, and the drop caused by the pdV work done as the 
outflow expands, respectively. The total luminosity. A, 
is the combined rate of energy deposition and dissipation 
within the outflow. 



A — LaGN — -^Comp — ill 



^diss : 



(13) 



where Lagn is the luminosity of the AGN contributing 
to the expansion of the outflow, Lcomp is the inverse 
Compton cooling off the CMB photons, Lne^ is the cool- 
ing due to two-body interactions. Lion is the cooling due 
to ionization, and Ldiss is the heat dissipated from col- 
lisions between the expanding shell and the IGM. We 
assume that the first two terms, AGN luminosity and 
inverse Compton cooling, dominates (following PMG07) 
and neglect the other terms. The Compton luminosity is 
given by 



^Comp 



2^ 

Is" 



arph 



fkT. 



^.jK-^) (l-cosf)(l+.)Vi^3, 

(14) 

(PMG07) where T^q is the temperature of the CMB at 
present. 

Magnetic fields thread AGN jets from sub-pc to kpc 
scales, but their detailed characteristics are not well 
known. Following FLOl, we make the simplifying as- 
sumption that during its activity period an AGN ejects 
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magnetic energy equal to a fixed fraction of the kinetic 
energy injected by the jets, Eb ~ sbEk- The fraction 
SB is taken as a free parameter, with a canonical value 
of es = 0.1. We also assume that the magnetic field (of 
strength B„i) is tangled, and exerts an isotropic pressure 
Pb = (1/3)^™/^''' affecting the dynamics of the out- 
flow. The relation between energy, pressure and volume 
is Eb = 3pgV. From this, we can derive the equation 
for the evolution of the magnetic pressure. 



Pb 



es^AGN 



47ri?3 [1 - cos (a/2)] 



R 



(15) 



where the last term comes from the conservation of mag- 
netic flux of a magnetic field frozen into the expand- 
ing outflow. At early times when the AGN is active, 
the first term in the right-hand side dominates, and 
we get Pb/Pt — After the AGN turns off, we 

set Lagn = (see below), and equation ([TS]) gives 
oc In our simulations we find that the thermal 

pressure always dominates over the magnetic pressure of 
the outflows. 

The combination of equations ([5))- p^ fully describe 
the evolution of the outflows. In Appendix A, we describe 
how these equations are solved in practice. 

In our simulations we allow each AGN to evolve 
through an active life when Zbir > z > ZoS- For 
this time period Tagn, we use the jet kinetic luminos- 
ity as the AGN luminosity in equations ([Tn|) and ([TS]) . 

^AGN = Lk = e^ibol- 

After the central engine has stopped activity (when 
z < Zos), it enters the post-AGN dormant phase. Then 
we set the AGN luminosity to zero (Lagn — 0) in equa- 
tions (|13p and (fT5|) . The gas inside the outflow is still 
overpressure d relative to the IGM , so the expansion con- 
tinues (e g.. iKronberg et all [200]] : iRevnolds et al.ll2002l 
PMG07, lBaraill2008f) . but the pressure drops faster since 
there is no energy input from the AGN. The outflow 
keeps expanding as long as its pressure prp + pg exceeds 
the external pressure p^ of the IGM. When Px+Pb — Px, 
the outflow has reached pressure equilibrium. After this 
point, the outflow simply evolves passively with the Hub- 
ble flow. 

2.6. Metal Enrichment 

The material transported into the IGM by AGN out- 
flows, for the most part, does not originate from the 
AGN itself, but from the interstellar medium of the host 
galaxy. Stellar evolution in the host galaxy have enriched 
the interstellar medium with heavy elements, or metals, 
and these metals will be carried into the surrounding 
IGM by the outflows. Using our simulations, we can cal- 
culate the distribution and concentration of metals in the 
IGM. Here we focus on the distribution of metals and the 
enriched volume fraction, that is, the fraction of the IGM 
by volume that has been enriched by AGN outflows. The 
concentration of metals in the IGM will be addressed in 
a forthcoming paper. 

To calculate the fraction of the total volume filled by 
outfiows, we cannot simply add up the final volumes of 
the outflows, and divide by the volume of the compu- 
tational box. Intergalactic gas enriched by outflows will 
move with time as structures grow, and therefore regions 
that were never hit by outflows might end up containing 



metals. We take this effect into account, by employing 
a dynamic particle enrichment scheme to quantify the 
cosmological enrichment history of the IGM. During the 
evolution of an outflow, we identify the particles that 
the outflow hits. These particles are flagged as been en- 
riched. This is done at every timestep for every outflow 
present in the simulation box at that time. The frac- 
tion of the volume of the box occupied by these enriched 
particles is then estimated, as follows. 

We divide the computational volume into = 256'^ 
cubic cells, and identify which cells contain matter that 
has been enriched by outflows. Notice that we cannot 
simply count the number of cells that contained parti- 
cles that have been flagged as enriched. This approach 
works in high-density regions, but fails in low-density re- 
gions where the cells are smaller than the local particle 
spacing. We solve this problem by using a Smoothed 
Particle Hydrodynamics technique. A smoothing length 
h is ascribed to each particle. We calculate h iteratively 
by requiring that each particle has between 60 and 100 
neighbors within a distance 1.7 h. We then treat each 
particle as an extended sphere of radius l.7h over which 
it is considered to be spread. The cells that are covered 
by one or more enriched particles are then considered 
enriched. This method works well both in low- and high- 
density regions. The total number of flUed cells, A^rich, 
gives the total volume of the box occupied by the en- 
riched particles. The fractional volume of the simulation 
box enriched by AGN outflows is then Ni-ich/N^. Notice 
that this is not done during the simulation itself. The 
simulation produces dumps at various redshifts contain- 
ing the positions and velocities of particles, as well as 
the flags which indicate which particles are enriched in 
metals. The calculation of the enriched volume fraction 
is then performed as post-processing. 

3. RESULTS AND DISCUSSION 

We performed a series 9 simulations by varying some 
parameters of our AGN evolution model. Table [1] sum- 
marizes the characteristics of each run. The flrst column 
gives the letter identifying the run. In columns 2, 3, and 
4, we list the opening angle of the outflows, the lifetime 
of the AGN's, and the kinetic fraction, respectively. The 
fifth column indicates if biasing was used when determin- 
ing the locations of the AGN's. The results are presented 
and discussed in the following subsections. For each run 
we plot the fraction of the cosmological volume enriched 
by the AGN outfiows. The last column of Table[T]lists the 
resulting metal-enriched volume fraction at the present 
epoch obtained for each run. 

3.1. Evolution of a Single Outflow 

Figure [3] shows the redshift evolution of a single 
anisotropic outfiow with opening angle a = 60°, born at 
z ~ 5.9. The different phases of the outfiow expansion 
are separated by vertical lines: active phase of the AGN 
on the left, post-AGN phase in the middle, and passive 
Hubble expansion on the right. The top panel shows 
the total luminosity (A = ^agn — ^Comp), the middle 
panel shows the comoving radius of the outfiow, and the 
bottom panel shows the pressures associated with the 
outfiow. 

Soon after birth, the outflow is driven by the active 
AGN, whose luminosity dominates over the Compton lu- 
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TABLE 1 

Parameters of the Simulation and Final Enriched Volume 
Fractions 



Run 


a ( ) 


Tagn (yr) 




Bias in Location 


AT / at3 / ~ r\\ 
Nnch/N^{z = 0) 


A 


180 


los 


0.10 


X 


0.80 


B 


120 


lO* 


0.10 


X 


0.82 


C 


60 


10« 


0.10 


X 


0.83 


D 


60 


10' 


0.10 


X 


1.00 


E 


60 


10^ 


0.10 


X 


0.75 


F 


60 


10« 


0.05 


X 


0.79 


G 


60 


10« 


0.01 


X 


0.71 


H 


60 


10« 


0.10 


V 


0.75 


I 


60 


10* 


0.10 


V 


0.65 



3.0-10" 




1 



Fig. 3. — Characteristic quantities for the evolution of a single 
outflow with lifetime Tagn = lO^yr, bolometric luminosity L^oi = 
3.3 X 10^-Lq, and opening angle a = 60°, as a function of redshift. 
Upper panel: total luminosity [eq. II13I I. A = Lagn ~ ^'Comph 
middle panel: comoving radius of outflow (-Rcomoving); lower panel: 
thermal pressure inside the outflow (pt, dashed), magnetic pressure 
(pb, dash- dot- dot- dot), total outflow pressure {pT + ps, solid), 
and pressure of the external IGM (px, dotted). The vertical dash- 
dotted lines in the lower two panels separate the different phases 
of expansion of the outflow: active, post-AGN overpressured, and 
the final passive Hubble evolution. 

minosity, and it grows rapidly in size. At the end of the 
active phase {z ~ 5.5), the AGN turns off and the only 
contribution to the total luminosity is the energy loss by 
Compton drag. The interior of the outflow is still over- 
pressured with respect to the external IGM by a factor 
{Pt ~^Pb)IPx ~ 400. So it continues to expand while 
its pressure falls faster (notice the change of slope in the 
bottom panel). Finally, when the total outflow pressure 
falls to the level of the external pressure, the expansion 
stops. From z ~ 1.8 the comoving radius of the outflow 
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o 
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> 0.4 

T3 
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0.0 



- a = 60° 
... a= 120° 
. - a = 180° 




Fig. 4. — Fractional volume {Nj-i^^^/Ng) of simulation box en- 
riched by AGN outflows as a function of redshift, for various open- 
ing angles: a = 60° {solid), 120°, {dotted), 180° {dashed). See ^?2\ 
for details. 

remains constant at 0.53/i~^Mpc in the passive Hubble 
flow phase. 

3.2. Opening Angle of Anisotropic Outflows 

To study the effect of varying opening angle of the out- 
flows on the enriched volume fraction, we performed runs 
A, B, and C with opening angles a = 180°, 120°, and 60°, 
respectively. Here a = 180° corresponds to the isotropic 
outflow case (the geometry used by FLOl, LG05, and 
other previous studies), and a — 60° represents the most 
anisotropic outflow we consider. 

Figure [4] shows the redshift evolution of the enriched 
volume fractions we obtained for different opening an- 
gles. The enriched volume fractions are small at z > 3, 
reach 0.1 at z ~ 2, and grow rapidly afterward to 0.4 at 
z ~ 1. At the present epoch, z = 0, 80% of the universe is 
enriched by AGN outflows. We do not find much depen- 
dence on the opening angle. More anisotropic outflows 
(smaller a) are found to enrich slightly larger volumes. 
Such a trend is counter-intuitive since, for a given ra- 
dius, an outflow with a = 180° occupies a larger volume 
than one with a = 60°. But according to our model 
prescription in §2.51 an outflow with a smaller a grows 
to a larger radius than a more isotropic one, because 
the energy is concentrated into a smaller volume, result- 
ing in a larger pressure. Still, isotropic outflows tend to 
have larger individual volumes than anisotropic ones, as 
PMG07 showed. 
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This effect is compensated by another effect: overlap- 
ping outflows. Because we locate AGNs in dense regions, 
they are strongly clustered, especially at low redshift. A 
dense region will eventually harbor many AGNs"* and 
their outflows will likely overlap. Isotropic outflows tend 
to overlap with one another more than anisotropic ones, 
for two reasons. First, for isotropic outflows the volume 
occupied by each outflow is larger, hence the probabil- 
ity of overlap is also larger. Second, anisotropic out- 
flows propagate along the direction of least resistance. 
As PMG07 showed, if several outflows originate from a 
common structure, like a cosmological filament or pan- 
cake, they tend to align themselves along the direction 
normal to the structure. This greatly reduces the amount 
of overlap, especially at small opening angle a. In our 
simulations with larger a, the outflows overlap more, re- 
sulting in a smaller enriched volume fraction. So the 
combined effect of smaller individual volumes and less 
overlap cause the more anisotropic outflows to enrich a 
slightly larger volume fraction of the IGM. 

Figures [S] and [H] show the evolution of the large-scale 
structure and the distribution of metals in a slice of co- 
moving thickness 2/i~^Mpc, for isotropic outflows (run 
A), and outflows with opening angle of 60° (run C). At 
z = 3.38, the metals are just starting to emerge from the 
dense regions that are hosting AGNs. As the simulation 
progresses, new AGNs are formed, producing additional 
outflows, while the ones already present keep expanding 
until they reach pressure equilibrium with the external 
IGM. By z = 1.05, a significant fraction of the volume 
is enriched, but metals are still concentrated near dense 
regions. By z = 0, the metals have spread into the low- 
density regions, and only the very-low-density regions, 
located far from any cosmological structure, have not 
been enriched. 

There is no striking difference between the two figures, 
though we can notice that voids are more enriched at 
z = 1.05 in run C than in run A. To illustrate this more 
clearly, we plot in Figure [7] the difference between the 
two enrichment maps. Red indicates regions that are 
enriched in run A but not in run C, while green indi- 
cates regions enriched in run C but not in run A. Since 
the enriched volume fraction is about the same for both 
runs, we expect comparable amounts of red and green, 
which is indeed the case. However, their distributions 
are significantly different. Looking at redshift z = 1.05 
(bottom-left panel of Fig. [7]), we find green areas all over 
the slice. Colored regions found inside deep voids are 
almost exclusively green, while red regions are predom- 
inantly located inside or near dense regions. Overall, 
anisotopic outflows enrich low-density regions more than 
isotropic outflows, at the expense of high density regions. 
This effect is seen most clearly at the intermediate red- 
shift z = 1.05. At earlier times, the enriched volume 
fraction is small in both simulations, making the differ- 
ence small, while at late times most regions are enriched 
in both runs, as the enriched volume fraction exceeds 
0.80. Figure [8] shows a zoom-in of the upper right re- 
gion, which illustrates the distribution of red and green 
areas relative to dense structures. 

3.2.1. Enrichment of Overdense vs. Underdense IGM 
* They might be active at different epochs, though. 



Fig. 5. — Shce (of comoving size 128 Mpc X 128 Mpc, 
and 2/i~^Mpc comoving thickness) off the computational volume 
for run A, showing the evolution of the distribution of metals. The 
areas shown in orange are enriched. The black dots represent the 
PM particles, and show the large-scale structures. 



Fig. 6. — Same as FigureO for run C. 



Fig. 7. — Differential enrichment map between runs A and C, 
calculated from Figures [5] and [6] Red areas show regions enriched 
in Run A but not in Run C; green areas show regions enriched in 
Run C but not in Run A. 



Difference A/C, z = 1.05, zoom-in 




Fig. 8. — Differential enrichment map between runs A and C, at 
redshift z = 1.05. This figure shows a zoom-in of the upper right 
region of the third panel of figure [71 

We investigate the relationship between metal enrich- 
ment by AGN outflows and the density of the regions 
being enriched, using runs A, B, and C. Figures [9l - fT3l 
show various density statistics of the enriched regions at 
five different redshifts, z = 3.38, 2.57, 2.0, 1.05 and 0, for 
3 opening angles, a = 60°, 120°, and 180°. All quantities 
are plotted as functions of the overdensity (ratio of the 
local density to the mean density, p/ Pmcan)- 

The top panels show the fractional number of cells in 
the whole grid (whether enriched or not) at a density 
p. The distribution resembles a skewed Gaussian, with 
more underdense volumes than overdense. The Gaussian 
becomes wider in time as large scale structures form in 
the cosmological volume, with further higher and lower 
densities being reached gradually. Most of the volume 
of the simulation box is underdense, and the peak of the 
distribution shifts to lower densities with time. 

The second panels show the fractional number of cells 
which are enriched by AGN outflows {N"'^^). The en- 
riched volume fraction is very small at early time, and 
increases with time, becoming significant from z ~ 2 — 3 
and afterward. At a given redshift, we see a slight 
trend of the peak of the resulting Gaussian distribution 
shifted to lower densities for increasing anisotropy {a = 
180°, 120°, 60°). This indicates that more anisotropic 
outflows enrich underdense regions of the IGM. This is 
most noticeable at redshift z = 1.05 (Fig. [T2|. 

The third panels show the ratio of the cell counts in 
the second to the top panel, i.e., the number of cells 
enriched as the fraction of the total number of cells. In 
high-density regions (p > Pmcan) more isotropic outflow 
enrich larger volumes than anisotropic ones, while in low- 
density regions (p < Pmcan) we flnd exactly the opposite 
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-3-2-10123 

log (P / Pmean) 

Fig. 9. — Density statistics at redshift z = 3.38 for runs A, B, 
and C. The quantities TVp, A^^''^'^, A^^pS and Af| = 256^ designate 
the number of cells (whether enriched or not) at a given density, 
the number of cells enriched by AGN outflows at a given density, 
the number of cells enriched by AGN outflows below a given den- 
sity threshold, and the total number of cells in the computational 
volume, respectively. Top panel: histogram of the density distri- 
bution of cells. Second panel: number of cells enriched by AGN 
outflows at a given density, divided by the total number of cells 
in the computational volume. Third: number of cells enriched by 
AGN outflows at a given density, divided by the total number of 
cells at that density. Bottom: number of cells enriched by AGN 
outflows below a certain density threshold, divided by the corre- 
sponding number for the isotropic case (ct = 180°). The linestyles 
indicate the different opening angles: a = 60° (dashed), 120° (dot- 
ted), 180° (solid). 

behavior. The effect is small, but most prominent at 
z = 1.05. At z = 0, dense regions are fully enriched for 
all values of a, while at z > 2 low-density regions are 
just starting to get enriched. 

The bottom panels are cumulative versions of the sec- 
ond panels, showing the number of enriched cells be- 
low a given density threshold {N"^^), divided by the 

corresponding number A'^"p'\8o isotropic outflows.^. 
In low-density regions [log(yo//9mcan) < 0], the ratio 
N<f/N<tiso strongly depends on the opening angle, 
at all redshift but especially at early time (2 > 2). The 

5 Note that Nlif" = N""^ in the limit of large p. Hence, the fact 
that the curves reach different asymptotes in the right hand side 
of the panels simply reflects the weak dependence of the overall 
enriched volume fraction on opening angle. 




z 1.0- 



-3-2-10123 

log (P / Pmean) 
Fig. 10. — Same as Figure[9] at redshift z = 2.57. 

most anisotropic outflow {a — 60°) enriches the largest 
underdense volumes, and the effect is stronger at lower 
densities. The effect is most prominent at earlier epochs 
(z 3.5 — 2), and decreases in intensity with time, but 
is still visible at the present epoch. We explain this 
trend as follows. At earlier times, when the AGN out- 
flows have just started to enrich the surrounding vol- 
umes, more anisotropic outflows enrich larger fractions 
of the lower density IGM, since these outflows expand 
along the direction of least resistance. As time goes on, 
large-scale structures form and grow, causing material 
to accrete from low- to high-density regions. This tends 
to even out the density of the enriched volumes, and as 
time goes on the effect of preferential enrichment of low- 
density regions by anisotropic outflows gets washed out. 
Such a trend in seen in our resulting density statistics, 
where the difference between the different opening angles 
is very small in Figure [13] at z = 0. 

3.2.2. Comparison with Fieri et al. 2007 

The bottom 3 panels of Figures [9HT31 are analogous to 
the plots in Figure 9 of PMG07. These authors consid- 
ered outflows driven by SNe explosions in dwarf galax- 
ies, and found that the anisotropy of the outflows had a 
major impact on the IGM enrichment in low-density re- 
gions. Outflows with opening angle a = 60° enriched the 
lowest-density regions of the IGM by a factor 3 — 4 more 



10 



0.15 
0.10 

0.05 

0.00 
0.08 

0.06 

0.04 

0.02 
0.00 

0.8 

Q. 

; 0.6 

0.4 

0.2 
0.0 

o 

CO 

? 2.0 
" 1.5 

Q. 

: 1.0 




z = 2 



■ 60° 
120° 
-180° 




-3-2-1012 

log (P / Pmean) 
Fig. 11.— Same as Figure l9lT0l at redshift z = 2.00. 
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Fig. 12.— Same as Figure [OTTTI at redshift z = 1.05. 



than isotropic outflows. In this study, where we consider 
outflows driven by AGN in active galaxies, we find that 
a similar trend of increasingly anisotropic outflows to en- 
rich lower-density regions, but the effect depends signif- 
icantly on the redshift. The trend we obtained is not 
as dramatic as in PMG07 at the present epoch (Fig. [13]) , 
but becomes more and more prominent at earlier epochs. 
We find that at z = 2.57 (Fig. [lOl), the factor by which 
outflows with a ~ 60° enrich low-density regions more 
than isotropic ones goes up to 1.7. At z = 3.38 (Fig. [9l), 
the factor goes as high as 2.3. 

Even though we have used the same prescription for 
the anisotropic outflows as PMG07, we point out the 
following differences between the two studies. PMG07 
did not perform an actual cosmological simulation. 
Instead, they generated an initial gaussian random 
fleld at high redshift, filtered that density field at 
various mass scales, and used the spherical collapse 
model to predict the collapse re dshift of each halo (see 
IScannapieco fc Broadhurstl200l1 and PMG07 for details). 
This semi-analytical model has the advantage of simplic- 
ity, but does not include an accurate treatment of halo 
mergers and accretion of matter onto halo, and further- 
more it undere stimates the level of clus tering of the halos 
(see, however. iPinsonneault et al.ll200 9). In this paper, 
we argue that the overlap between outfiows, which re- 
sults from the clustering of halos, is more important for 



TABLE 2 

Runs with different active lifetimes Taqn 



Run 


Tagn (yr) 


Population Size 


A^rich/A^l {z = 0) 


D 


10^ 


15 279 029 


1.00 


C 


10« 


1 535 362 


0.83 


E 


109 


162 228 


0.75 



isotropic outfiows than anisotropic ones, and that ac- 
cretion of low-density, metal-enriched matter onto dense 
structures tends to erase the effect of anisotropy. As a 
consequence, we find that the effect of anisotropy is still 
important at high redshifts, but greatly reduced at low 
redshifts. 

3.3. AGN Lifetime 

The simulations A, B, and C presented in the previ- 
ous section all assumed that AGN activity lasts Tagn = 
lO^yr. To study the effect of varying Tagn, we per- 
formed runs D and E, using AGN activity lifetimes of 
2agn = lO^yr and lO^yr, respectively. These values cor- 
respond to the lower limit and upper limit on Tagn used 
by LG05. Simulations C, D, and E use the same opening 
angles, a — 60°, and the same kinetic fraction, e^^ — 0.1, 
and differ only in the value of Tagn ■ We neglect any re- 
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of few xlO^yr is considered as the typi- 
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Fig. 13.— Same as Figure l9lT2l at redshift z = 0.00. 
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Fig. 14. — Fractional volume (AfjichZ-'Vg) of simulation box 
enriched by AGN outflows as a function of redshift, for different 
AGN activity lifetimes: Taqn = 10 ''yr {dashed), lO^yr (solid), 
lO^yr {dotted). See jOl for details. 

peated activity (which might happen because of duty cy- 
cle) of the same AGN after it has become inactive, as did 
LG05. Considering duty cycles would have a small effect 
in our results, because of the way we model the outflow 
expansion and evolution (i i2.5[ our outflows continue to 
remain in the simulation box even after the central AGN 



cal activity lifetime of the SMBH a t the cen ters of 
active galaxies in other studies (e.g., lYu &: Lul [20081 . 
A lifetime of 5 x lO^yrs has b een used in t heoreti- 
cal modeling o f radio galaxies bv lBlundell et al.l (|l999); 
iBarai fc Wiital (l200'<f). Ob servations of X-ray activity in 
AGN (iBarger et al.ll2001h SPSS optical studies of ac- 
tive galaxies (jMiller et al.l 120031). and black ho le demo- 
graphics arguments (e.g.. [Marconi et al.ll2004 ) all sup- 
port an AGN activity li fetime of ~ lO^yr or more (also, 
iMcLure &: Dunlopll2004 l. 

Table [2| gives the parameters and results of runs C, D 
and E. The first and second columns give the run name 
and the value of Tagn used, respectively. Column 3 gives 
the size of the AGN population, that is, the total number 
of sources obtained from the QLF f iJ2.2p in our cosmo- 
logical volume over the Hubble time. The last column 
lists the resulting metal-enriched volume fractions at the 
present epoch. There are ^ 10 times more AGNs gener- 
ated for a decrease of Tagn by a factor of 10, since the 
QLF has to remain constant. This change in the num- 
ber of sources makes a large difference in the resulting 
enriched volume fractions. 

The redshift evolution of the enriched volume fractions 
for different active lifetimes are plotted in Figure [T4| The 
largest AGN population (Tagn = lO^yr) enriches signif- 
icantly larger volume than the other two populations. In 
this case (run D) the fractional volume starts to become 
appreciable from z ^ 3.5, and grows afterward, such that 
100% of the volume of the box is enriched by z ~ 0.7. 
The other two populations, with Tagn — 10^ and lO^yr, 
enrich 0.83 and 0.75 of the volume at z = 0, respec- 
tively. This result, that AGNs with shorter lifetimes en- 
rich larger volumes than those with longer lifetimes, is 
opposite to what LG05 and Barai (2008) obtained. We 
explain this in the following, stressing the point that the 
method we employed to enrich our cosmological volumes 
is quite different from what other studies have used. 

We use a dynamic particle enrichment scheme (i i2.6p . 
whereby we marked particles lying inside outflow vol- 
umes as enriched, assign a smoothing length to each of 
them, and calculate the volume of the box being filled 
by these en riched particles. In other studies (LG05 and 
lBarai|[2008[ ) the enriched volume fraction is calculated by 
finding the fraction of the computational volume which 
is located inside the geometrical boundaries of the out- 
flows. In our approach, particles that are metal-enriched 
by outflows are allowed to move, and might move to re- 
gions not intercepted by an outflow. This increases the 
chance of enriching a larger fraction of the total volume. 
The trend heightens when the number of outflows in- 
crease; as more AGN are born they are placed in distinct 
locations distributed throughout the box (according to 
§2.3p . and they could potentially enrich larger volumes. 
This causes the enriched volume fraction to increase with 
increasing population size, even though the lifetime of 
each outflow is shorter. 

3.4. Kinetic Fraction 

In runs F and G, we consider different values of the ki- 
netic fraction (the fraction of the AGN bolometric lumi- 
nosity which converts to the kinetic luminosity of the jets, 
H2.5L ck = 0.05 and 0.01, in addition to our standard 
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Fig. 15. — Fractional volume (Nj-i^i^/N^) of simulation box 
enriched by AGN outflows as a function of redshift, for different 
kine tic f ractions: ex = 0-1 (solid), 0.05 (dashed), 0.01 (dotted). 
See jO for details. 

value = 0.1 (run C). Comparing with other studies, 
FLOl used a value e^j- ^ 0.1.^ LG05 first experimented 
with a value = 0.1, and found an enriched volume 
fraction that was too large to agree with observations, so 
later they chose = 0.01 as their fiducial value. Hence, 
we are sampling the range of values that have been con- 
sidered in other studies. 

The resulting enriched volume fractions are plotted in 
Figure [T51 as a function of redshift, for the different ki- 
netic fractions. We find the same general trend as in 
runs A - E, namely that the enriched volume fraction is 
small at z > 2, and then grows rapidly afterward up to 
the present. At z = 0, fractions of 0.83, 0.79, and 0.71 of 
the volume are enriched by outflows with = 0.1, 0.05, 
and 0.01, respectively. 

Using their model for the AGN outflows, LG05 found 
that the cosmological volume is completely filled (100%) 
by z ~ 2 when using ck — 0.1 or 0.05, and with 
eK — 0.01, their volume filling fraction at z = is 80%. 
The enriched volume fractions we obtained with = 0.1 
and 0.05 are comparable to those obtained by LG05 with 
€k = 0.01. With €k = 0.01, we obtained an enriched 
volume fraction ~ 10% smaller than that of LG05. We 
find that decreasing the kinetic fraction by a factor of 10 
reduces the volume enriched by ~ 12%. This is compa- 
rable to a corresponding « 20% reduction in the volume 
filling fraction found by LG05 (their Fig. 3), when ex is 
reduced by a factor of 10. It may seem surprising that 
increasing by a factor of 10 the energy driving the ex- 
pansion of the outfiows has only a ~ 14% effect on the 
final results. With a larger kinetic fraction, outfiows are 
certainly larger, but also overlap more with one another. 
The lowest-density regions, located far from any AGN, 
still manage to escape enrichment. 

3.5. Bias in AGN Location 

AGNs are observed to be clustered in the universe, 
with the clustering aniplitude increasing with red- 
shift (jPorciani et al.ll2004l : ICroom et al.l[2005l : IShen et all 
[200l . We study the effect of bias in the distribution of 

FLOl quote a value = 1.0, but they expressed the quasar 
kinetic luminosity in terms of the B-band luminosity instead of the 
bolometric luminosity. 
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Fig. 16. — Fractional volume (Nj-i^i^/N^) of simulation box 
enriched by AGN outfiows as a function of redshift, The linestyles 
indicate different methods of AGN location: Default run C (solid), 
run H with bias parameter ip = 2, as in LG05 (dotted), run I, biased 
to densest peaks (dashed). See i]3.5l for details. 

AGNs by considering two different methods for spatially 
locating the AGNs inside the computational volume, in 
addition to our original method ( §2.3p used in runs A-G. 

In run H, we use a location biasing condition similar 
to that used by LG05, but restricted to cells containing 
a local density peak (this is required by our method for 
finding the direction of last resistance as discussed in 
H2.3\i . We calculate the probability of each peak cell of 
hosting an AGN as 

'P=^v^ . (16) 

where p is the filtered density (see iJ2.3p of the cell, A^pcak 
is the number of available density peaks (not containing 
an AGN already) inside the computational volume, and 
tp is the bias parameter, which is denoted by a in equa- 
tion (2) of LG05. The summation in the denominator 
of equation (flB]) is over all the density peaks at the rel- 
evant timestep. We use a constant bias value 7/1 = 2, 
as did LG05 in their fiducial run. At each timestep, we 
use a Monte Carlo rejection method to locate the AGNs 
randomly, with a probability given by equation 

In run I, we bias the AGNs in favor of the highest- 
density peaks inside the cosmological volume. At each 
timestep, we find all the peak cells fi )2.3p . and then sort 
them in descending order of their density. The new AGNs 
(those born in that timestep) are located in the peak 
cells which have the highest densities. We start with 
the highest-density peak and locate one AGN (selected 
randomly from the ones just born) in it, go to the next 
highest-density peak cell and repeat the process, until 
all AGNs have been assigned locations. Notice that the 
methods used in simulations C and I effectively corre- 
spond to bias parameters "0 = and ^ 00, respec- 
tively. 

Figure [16] shows the redshift evolution of the enriched 
volume fraction for the different location methods. The 
enriched volume fraction decreases with increasing bias, 
at all redshifts. At z = 0, the enriched volume fraction 
reaches 0.83 for the case without bias {ip — 0, run C), 0.75 
for the intermediate case {tp — 2, run H), and 0.65 for the 
maximum bias case {tp 00, run I). Such a behavior is 
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expected, as with the introduction of gradual bias AGNs 
are increasingly clustered in high-density regions, so their 
outflows tend to overlap. Biasing to the highest density 
peaks reduces the present enriched volume fraction by 
18% compared with the default case. 

4. SUMMARY AND CONCLUSION 

We have implemented a semi-analytical model of 
anisotropic AGN outflows in cosmological N-body sim- 
ulations. The AGNs are placed at the location of local 
density peaks, and the outflows are allowed to expand 
anisotropically from those peaks, along the direction of 
least resistance, with a biconical geometry. Each AGN 
produces an outflow from its birth, which at first ex- 
pands rapidly in the active- AGN phase for a time Tagn 
until the AGN turns off, and then continues to expand 
slowly due to its overpressure. Finally after coming to 
pressure equilibrium with the external IGM the outflow 
simply follows the Hubble flow. The outflows carry with 
them metals produced by stars in the host galaxy, and 
enrich the regions of the IGM that they intercept. We 
used this algorithm to simulate the evolution of the large- 
scale structure and the propagation of outflows inside a 
cosmological volume of size (128 Mpc)^, from initial 
redshift z = 24 to flnal redshift ^; = 0, in a concordance 
ACDM model. We performed a total of 9 simulations, 
varying the opening angle of the outflows, the active life- 
time of the AGNs, the kinetic fraction, and the method 
used for locating the AGNs inside the cosmological vol- 
ume. 

(1) The enriched volume fractions (fraction by volume 
of the IGM enriched in metals by the outflows) are small 
at z > 2.5, and then grow rapidly afterward up to z = 0. 
In our simulations with different parameter values, we 
found that AGN outflows enrich 0.65 — 1.0 of the total 
volume of the universe by the present epoch. 

(2) The enriched volume fractions do not depend sig- 
nificantly on the opening angle of the outflows. More 
anisotropic outflows (smaller a) are found to enrich 
slightly larger volumes than more isotropic ones, because 
the anisotropic ones grow bigger in radius and have less 
overlap. Increasingly anisotropic AGN outflows enrich 
lower-density volumes of the IGM, to an extent that de- 
pends on redshift. The trend is most prominent at ear- 
lier epochs, and not very dramatic but still visible at 
the present. Outflows with opening angles a = 60° en- 
rich the underdense IGM more than isotropic outflows 
by a factor of 2.3 sX z — 3.38, and 1.7 at z = 2.57. 
Initially, high-density regions get enriched by outflows, 
since the AGNs producing these outflows are located pre- 
dominantly in these regions. Eventually, expanding out- 



flows reach low-density regions of the IGM, and these 
regions also get enriched. This effect is larger for more 
anisotropic outflows, since these outflows expand along 
the directions of least resistance, reaching the low-density 
regions sooner and more easily. That enriched matter, 
located in underdense regions at earlier epochs, gravita- 
tionally accretes into higher-density regions with time, 
as large-scale structures grow. As a result, the effect 
of preferential enrichment of underdense IGM by more 
anisotropic outflows gets washed out with time, such that 
the difference is quite small at the present, even though 
it is quite significant at redshifts z > 1. 

(3) Reducing the active lifetime Tagn of the AGNs re- 
sults in larger enriched volume fractions. Any reduction 
in Tagn must be accompanied by a corresponding in- 
crease in the number of AGNs, such that the observed 
QLF remains the same. When the number of sources 
increases, these sources tend to be distributed more uni- 
formly throughout the computational volume, resulting 
in a more efficient enrichment. With Tagn = lO^yr, 
100% of the volume is enriched by redshift z = 0.7. 

(4) The enriched volume fractions for kinetic fractions 
= 0.1,0.05 and 0.01 are 0.83,0.79 and 0.71, respec- 
tively, at the present epoch. This is consistent with the 
results reported by LG05. Increasing e^^ results in larger 
outflows, but the main consequence is to increase the 
level of overlap, rather than the enriched volume frac- 
tion. 

(5) We varied the prescription for locating AGNs in the 
computational volume by introducing a bias parameter 
which favors high-density peaks. This increases the level 
of clustering of the AGNs, and when AGNs are more clus- 
tered, there is more overlap between outflows, resulting 
in a lower enriched volume fraction at all redshifts. The 
maximum effect is an 18% decrease at z = when going 
from an unbiased distribution to a maximally biased one. 

This paper focused on the metal-enriched volume frac- 
tion by the outflows. In a forthcoming paper, we will fo- 
cus on the metal abundances in the IGM resulting from 
anisotropic AGN outflows, and their observational con- 
sequences. 
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APPENDIX 



NUMERICAL SOLUTION FOR THE OUTFLOW 



In this appendix, we describe the technique used for solving the equations describing the evolution of the outflows. 
The presentation closely follows the one used in the appendix of PMG07. 
The equations governing the evolution of the outflow are 
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with the initial conditions i? = at < = 0. We used equations ([5]) and pT|) to ehminate Md and px in equations (|A1[) 
and (|A2p . Notice that sinc e A/^ cx 1 — cos(a/2), the dependencies of the opening angle a cancel out in the first and 
last terms of equation (jAl[) . The only remaining term in that equation that depends on a is the term —GMg/2R^ . 
The most important dependencies on a are found in equations (jA3[) and (|A4p . The energies (thermal and magnetic) 
are injected into smaller volumes, leading to larger pressures. 

Many terms in equations (|Aip ~ (jA4p diverge at < = 0, making it impossible to find a numerical solution in this form 
(these divergences occur because we neglect the finite size of the source; in the real universe, many terms are large, but 
finite, at a radius R that it small, but nonzero). To improve the behavior of the equations at early time, we perform 
a change of variable. To find the appropriate change of variable, we take the limit i — > in equations (|Aip - (jA4p . and 
keep only the leading terms. We get 
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where the subscripts i indicates initial values a.t t — 0. Notice that in equation (jASp . we used Lagn 
t ^ 0. We can easily show that the solutions of equations (jA5p - (IA8P are power laws. 
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In the limit < ^ 0, the six functions S, Ns, qrj^, Qg, U, and Os vary linearly with t. We now eliminate the functions R, 
Ms, prp and in our original equations (|Aip - (|A4p using equations (IA14p - ()A19p and get 

145" ATrS^iq^ + qB-qx) f. a\ GNL^it^^'"" nu'^St 



U-- 



91/ 



t 5t 



(A20) 



15 



<1t = 



4jV, 
5t ' 

2G 
At3 



(l-cos|) 



U 

t3 



5" = 



27r53[l -cos(a/2)] 



bqrpU 
St 



25 _ 
5t 



l2~ 



5i 



HSt 



47r53[l -cos(a/2)] 
[/ 
T' 

where =p^t^/^ in equation (jA20p . 

These equations are completely equivalent to our original equations (jAl[) - (|A4[) . but all divergences at i = 
been eliminated. These equations can therefore be integrated numerically using a standard Runge-Kutta alg 
with the initial conditions U = S = Ng = qj, = q^ = Q &t t = Q. However, before doing so, it is preferable to 
the equations in dimensionless form. We define 



a = 

Ml 



I Hot , 
HqS 

' C ' 
HqU 

' C ' 

_ GNs 

_ GOs 

GqT 
'C^o 

, GqB 
, Gq^ 

H 

A 

'aI' 



(A21) 

(A22) 
(A23) 
(A24) 

have 
orithm, 
rewrite 

(A25) 
(A26) 

(A27) 
(A28) 
(A29) 
(A30) 
(A31) 
(A32) 
(A33) 
(A34) 
(A35) 
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where Os = ridNs/dr) in equation (jA36|) . In the Umit r — > 0, where many terms take the form 0/0, the derivatives 
reduce to dS/dr = dtj/dr = 1, dg^/dr = 2mb,,/|r .y2007r(l + 65/2), dq^/dr = 2117t, J| ./2007r(l + 2/€b), and 

diVs/dr = nb,ipH,i (1 - cos a/2) /2. 

The quantity /a appearing in equations (|A38|) and (|A39|) depends on the luminosity Lp omp, w hich is given by 
equation (Hi]). We chminate Pt a nd R in equation JTl]), using equations (|A14|) . (|A16p . (|A25p . (|A26p . and (|A30p . then 
ehminate C using equation (|A13p . We get 
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